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^ ! Abstract 

A code that impiements Einstein equations in the characteristic formulation 
^ in 3D has been developed and thoroughly tested for the vacuum case. Here, 

we describe how to incorporate matter, in the form of a perfect fluid, into 
the code. The extended code has been written and validated in a number of 
\ cases. It is stable and capable of contributing towards an understanding of a 

number of problems in black hole astrophysics. 
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I. INTRODUCTION 



A code based on the characteristic formulation of numerical relativity has been developed 
for the general 3D problem [Q . The code computes the gravitational field of the full Einstein 
equations for a vacuum spacetime between some inner timelike worldtube T and future null 
infinity. The code has been tested with the inner worldtube T being the (past) event horizon 
of a Schwarzschild black hole, with incoming gravitational radiation scattered off the black 
hole. The tests have included both the linear and highly nonlinear regimes, and show the 
robustness of the characteristic formulation in the study of radiative problems. Here we 
consider whether there are real astrophysical problems to which the characteristic code 
alone could be applied. 

Astrophysical problems that can be tackled by the code require the existence of a natural 
inner worldtube T on which boundary data is known. This would include the (past) event 
horizon of a black hole. In principle the black hole could be of Kerr-Newman type, but 
at present suitable boundary and initial data (for the characteristic formulation) are known 
only in the Schwarzschild case, and also in some multi-hole vacuum spacetimes || . However, 
if matter is incorporated into the code, then the range of possible astrophysical applications 
can be greatly extended. It would then in principle be possible to compute, in full nonlinear 
general relativity, the gravitational field and matter flow of a black hole accreting dust and 
gas, and of a black hole capturing polytropic or more realistic approximations to neutron 
stars. This paper does not solve any of these important astrophysical problems, but rather 
it lays the groundwork for doing so. 

The incorporation of matter into a characteristic code has been discussed previously. An 
axisymmetric code with matter has been used to consider a problem in cosmology || ; and a 
spherically symmetric matter code, that includes Cauchy-characteristic matching, has been 
reported recently In contrast, the present work has no symmetry requirement. 

The paper restricts attention to matter in the form of a perfect fluid, by which we 
mean that the pressure p is a function only of the density p. In section [IT] we summarize 
known results, and introduce our notation, for the characteristic formalism and perfect fluid 
evolution. Next, in section |TJ we derive the fluid evolution equations in our characteristic 
coordinates; we also find the additional terms (compared to the vacuum case) that appear in 
the Einstein equations. Section [TV] describes details of the numerical implementation of the 
fluid evolution equations (the additional terms in the Einstein equations are straightforward 
to code). The resulting matter plus gravity code has been run on a variety of test problems, 
and the results are described in section |V|. In all cases, the code is found to be stable 
and convergent. We end with a Conclusion (section |V7j), and an Appendix that contains 
expressions from section |TTJ that are rather long. 



II. FORMALISM 

A. The null cone - previous results 

The formalism for the numerical evolution of Einstein's equations, in null cone coordi- 
nates, is well known (see also J7| P~D|] ) . Nevertheless, for the sake of completeness, we 
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give here a summary of the formalism, including the necessary equations. We will use the 
notation and language of [0. 

We use coordinates based upon a family of outgoing null hypersurfaces. We let u label 
these hypersurfaces, x A (A = 2,3), label the null rays and r be a surface area coordinate. 
In the resulting x a = (u, r, x A ) coordinates, the metric takes the Bondi-Sachs form p~0|JTT 



ds 2 = - ^e 2/3 (l + y) - r 2 h AB U A U B ^j du 2 - 2e 2(5 dudr - 2r 2 h AB U B dudx A 

+ r 2 h AB dx A dx B , (1) 

where W is related to the more usual Bondi-Sachs variable V by V = r + W; and where 
h AB h B c = $c an d det{h AB ) = det(q A B), with q AB a unit sphere metric. In analyzing the 
Einstein equations, we also use the intermediate variable 

Qa = r 2 e- 2 ?h AB U B . (2) 

We work in stereographic coordinates x A = (q, p) for which the unit sphere metric is 

q AB dx A dx B = -^{dq 2 + dp 2 ), (3) 

where 

P = l + q 2 +p 2 . (4) 
We also introduce a complex dyad q A defined by 

ff*=f(M) (5) 

with i = y/—l. For an arbitrary Bondi-Sachs metric, h AB can then be represented by its 
dyad component 

J = h AB q A q B /2, (6) 

with the spherically symmetric case characterized by J = 0. The full nonlinear h AB is 
uniquely determined by J, since the determinant condition implies that the remaining dyad 
component 

K = h AB q A q B /2 (7) 

satisfies 1 = K 2 — J J. We also introduce spin-weighted fields 

U = U A q A , Q = Q A q A , (8) 

as well as the (complex differential) eth operators 9 and 8 (see [0 for full details). 

The Einstein equations G a b = decompose into hypersurface equations, evolution equa- 
tions and conservation laws. Naturally, the equations will require additional terms to allow 
for the presence of matter, and the extended equations are given later in section |U1B|. Here 
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we just note that the hypersurface equations form a hierarchical set for /3 r , (r 2 Q) jr , U >r and 
Wy, and the evolution equation is an expression for (rJ) tUr . 

The remaining independent equations are the conservation conditions, but they will not 
be needed here. 

The null cone problem is normally formulated in the region of spacetime between a 
timelike or null worldtube T and Z + , with initial data J given on the null cone u = in this 
domain. Boundary data for j3, Q, U, W and J is also required on T. The metric variables 
used remain regular at Z + , and we represent X + on a finite grid by using a compactified 
radial coordinate x = r/(l + r). 



B. Perfect fluid - previous results 



The description of a perfect fluid is well understood (e.g. [|13|,|14]]). The stress-energy 
tensor is 



- ab 



(p + p)v a v b + P9ab (9) 



where p is the density, p is the pressure, v a is the 4-velocity and g ab is the metric. In the 
cases considered here, the pressure depends only on the density, i.e. 

p = p(p). (10) 

The matter evolution equations follow from the conservation law T ab - b = 0, and are 

P,aV a + (P + P)V% = 0, ( P + P )v a;b V b + (S b a + V a V h )p, b = 0, (11) 

which may be rewritten as 

P, a v b g ab + (p + p)(v a>b - T c ab v c )g ab = (12) 
M a = {p + p)(v a>b - r c ab v c )v d g bd + p A + v a v cPib g bc = 0. (13) 

The numerical implementation of the fluid equations coupled to G.R. has been investigated 
primarily within the 3+1 or Cauchy framework. In this formulation, the spacetime is foliated 
by a sequence of spacelike surfaces, initial data is given on an initial surface and the evolution 
equations are used to compute the future. The numerical investigation of the "G.R.-Hydro" 



problem started in the early 70's by Wilson [15| and since then it has received considerable 



attention. The difficulty of modeling these equations has spurred the development of sophis- 
ticated techniques to deal with the diverse idiosyncracies of the problem. Thus, artificial 
viscosity techniques, total variation diminishing flux limiters, shock-capturing schemes, etc. 
are actively employed to aid in the numerical modeling; refer to |IB| for a recent review. 

Another problem one faces when attempting a numerical simulation in the 3+1 formalism 
is that an artificial "outer boundary" has to be included at some radius in order to deal 
with a finite grid. This introduces spurious reflections that spoil long term evolutions. In 
the characteristic formulation, on the other hand, one can use Penrose's compactification 
techniques to include infinity in the numerical grid. Additionally, being able to access null 
infinity allows one to obtain physical quantities, like the radiation given off by the system 
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and the total mass, unambiguously. There is not as much experience with the G.R.-Hydro 
problem in the characteristic formulation, but the results obtained so far are encouraging, 
indicating the characteristic formulation of G.R. might be a valuable tool to study a variety 
of astrophysical problems. 

III. PERFECT FLUID IN THE NULL CONE FORMALISM 

A. Perfect fluid equations 

In order to proceed further we represent the angular part of the 4-velocity by means of 
the complex quantity 

p 

V an g = q A VA = — (>2 + IV3), (14) 

where the ang suffix is introduced to avoid confusion with the Bondi-Sachs metric variable 
V. The matter evolution equations are then: 

• equation (|T^); 

• equation ( |13|) in the form 

M x = 0, q A M A = 0. (15) 



We also introduce the notation 



and then write 



1 dp na , 

Pp = ^—3- ( 16 ) 
p + pdp 



P,a= Pp(P + P)p,a- ( 17 ) 



The result of doing this is that, in the matter evolution equations, there is no explicit division 
by (p + p) (which could be zero). From a numerical point of view, it is possible to write the 
procedure for computing p p so as to ensure its good behavior in the low density limit. 

The matter evolution equations have been calculated using Maple; they are rather long 
and are given in the Appendix. The denominators are important in determining whether 
there may be singular behavior, and here we note that the form of the equations is 

1 

P ' U = rW(J-l) P (18) 



vi,u = 3 (dp F, (19) 

r3 Mf P - 1) 
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(Iliq.U o 

r z v± 



f v . 



(20) 



Formally, the equations ([18]) and (|H5|) could be singular if were to be 1; but this would 
correspond physically to the velocity of sound being equal to that of light. We do not 
concerns ourselves with those cases at the moment and defer its treatment for a later work. 
Also, the equations (0), ( |19| ) and (|2~0|) could be singular if v\ were to be zero; yet, starting 



from —1 



V a Vb9 



ab 



Vie 



-WfYi Vl -2v -2U A v A )<-l, 



(21) 



which can never be satisfied if V\ = 0. Thus, from an analytic point of view, the equations 
([[ID, (0) and ( p0|) form a well-behaved set of evolution equations. 

Finally, again using the condition — 1 = g ab v a Vb, we obtain an expression for the remaining 
velocity component Vq 



v 



e^(2KV ang V ang - JV a j g - JV a j g + 2r 2 ) + 2rv l {Vv l - rUV ang - rUV ang ) 

Avir 2 



(22) 



B. The Einstein equations 

The introduction of matter also amends the Einstein equations. We write the equations 

as 

Rob = 87r(T ab - ^g ab T), (23) 

and note that 

T = -{p + p) + 4p = 3p - p. (24) 
Thus the Einstein equations for a perfect fluid are 

3^) — p 

Rab = 8?r(p + p)v a v b + g ab (p — )) 

P — P 

= 87r((p + p)v a v b + g a b—^— )• (25) 

In the expressions below, Np, Njj, Nq, Nyy and Nj represent the nonlinear aspherical 
terms (in a sense specified in ||). These quantities are quite long and are not repeated here 
since they have been given explicitly in Using Maple we have confirmed that 



-Rn in equation (p5|) leads to 

# r = Np + 2irr(p + p)(v 1 ) 2 . (26) 
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RiAq A in equation ( p5|) leads to 

(r 2 Q), r = -r 2 (Bj + QK), r + 2r 4 5 {r~ 2 (3) r + N Q + 16vrr 2 (p + pHKng- (27) 

The equation for {7 is a definition so the presence of matter does not change it: 

p 2/3 

U r = —(KQ-JQ). (28) 

Rab^ AB in equation fl2~5|) leads to 

W, r = ^e 2/3 K - 1 - e p me? + V 2 (r 4 (3C7 + 3*7) ) r + N w 
-Ane^((p + p)(Ja^iU - \{JVl g + JKn S )) + (P " pV\ (29) 



where 



n = 2K- + i(3 2 j + s 2 j) + -^(ajsj - gjsj). (30) 



RabQ A Q B in equation fl25|) leads to 

2 M). ttr - (r^V (*"■/),,.) = -r" 1 (r 2 5£/) r + 2r~ 1 e /3 5V - (r" 1 ^ r J + Nj 
+—(V 2 ng (p + P )+r 2 J(p-p)). (31) 



C. Summary 

The data required on the initial null cone is: 

J P Vi V ang (32) 

and these constitute the set of evolution variables. The auxiliary variables may then be 
determined on the initial null cone, and they are found in the following order: p from the 
equation of state, (3 from equation (p6|), Q from equation (pTj), U from equation (p8|), W 



from equation fl29|) , and Vq from equation fl22|). The evolution equations (|3l|), (|18D , (|19|) and 



(BO) may now be used to find J, p, Vi and (in that order) on the "next" null cone. 

In order to have a properly specified problem, we will also need boundary data on the 
timelike worldtube 17 For the gravitational variables this data is (3, Q, U, W and J, and 
for the matter variables we will need p, v± and V ang . 
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IV. NUMERICAL IMPLEMENTATION 



We constructed a set of algorithms to solve equations (E(MJI]). In discretizing the equa- 



tions we follow closely the strategy developed for the vacuum case |y. We introduce a 
compactified radial coordinate x = r/(l + r) and define the numerical grid with coordi- 
nates (u n ,Xi,qj,pk) = (nAu, 1/2 + (i — I) Ax, —1 + (j — 3)Aq, —1 + (j — 3)Ap) (where 
2 Ax = l/(N x — 1), and Aq = Ap = 2/(N% — 1) ) . Using finite differences to discretize the 
equations, we center the derivatives at (n + 1/2, i — 1/2, j, k). 

The evolution equation is treated as in [|lj where the right hand side is modified to include 
the matter terms. Thus, the matter variables at (n + 1/2, i — 1/2) are evaluated by taking 
the average between the values at (n+l,i — 1) and (n,i), while the radial derivatives are 
obtained from the average of the corresponding values at (n + 1, i — 3/2) and (n, i + 1/2). 

Next, we proceed to integrate the evolution equation of the matter variables by a simple 
iterative method which results in a second order in space, second order in time scheme, as 
follows. First, note that the matter equations can be schematically written as 

9,u = F, (33) 

then, a direct second order discretization is obtained by 

g^ 1 = 9? + AuFr 112 ■ (34) 

Note that F depends on the matter fields, thus without having at hand the value of g™ + , 
F can not be directly evaluated at the midpoint. However, the value of _p™ +1//2 can indeed 
be obtained by a straightforward iteration. In the first step we set _F™ +1 / 2 = FT 1 an d use it 
to evaluate the right hand side of (|34D ; then, the obtained approximation to is used to 
obtain a better approximation to _F. n+1//2 an d so on. This procedure is repeated a sufficient 
number of times to ensure second order convergence of the algorithm. However, care must 
be taken at the horizon where fields diverge which will spoil the numerical modeling of the 
problem. If we had at hand a solution that we could use to evaluate the fields near the 
horizon, we could just integrate away from it. Unfortunately, this is not the case, and a 
different strategy must be employed. 

In this work, where we mainly investigate the suitability of the characteristic formulation 
as a mean to simulate matter coupled to G.R.; we resort to evolving the matter equations 
from the first point outside the worldtube to null infinity. Hence, the previously described 
algorithm, where F is evaluated at the i-th point, can not be used since F involves radial 
derivatives. Hence, for the first point we use a different algorithm by employing a one sided 
scheme taking backward differences to evaluate F. This scheme is only first order accurate 
and its resulting discretization error does not match smoothly to the one obtained with 
the second order scheme; thus, we do not expect to obtain global second order accuracy in 
our numerical integration. Several other alternatives are being investigated to obtain global 
second order accuracy but we defer their implementation to a future work. 

Finally, the hypersurface equations are discretized as in PJ, where the right hand sides 
now include the matter variables evaluated at (n + 1, i — 1/2) (which is straightforward 
having the matter fields values at (n + obtained in the previous step). 
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V. TESTS AND RESULTS 



We test the code by considering an initial localized distribution of matter around a 

Schwarzschild black hole with mass taken to be M = 1. The gravitational initial data is 
taken as 

J(u = 0,r,x A ) = 0. (35) 

For a non-spherical initial distribution of matter, this condition is in a sense unphysical in 
that it will introduce spurious incoming gravitational radiation [17|,18]; nevertheless it is 
simple and is suitable for code testing. The gravitational boundary data on the worldtube 
T, which we take to be the (past) event horizon of the black hole at r = 2, is 

P = 0, Q = 0, U = 0, W=-2, J = 0. (36) 

The initial data for the tests include two different distributions of the matter 

• Spherical shell that falls radially into the black hole, 

• Localized blob of matter falling radially towards the black hole. 
The tests are performed for two different equations of state 

• Dust (p = 0) 

• Fluid with p oc p 1A . 



A. Initial and boundary data for the matter 

We assume 

,>(,/ = 0.r..r l ) = { V 2<r-R*)J V 2(R b -r)J (3?) 

otherwise, 

If G(x A ) = 1, p describes a spherical shell of matter between r = R a and r = Rf,, and 
centered about r = 0. We also consider the case with G defined as a localized gaussian-type 
function 

f (q 2 +p 2 - aO 4 if ? 2 + p 2 < v 

G(x A ) = I (38) 
[ otherwise, 

which we use to describe an initial "blob" of matter. To provide initial data for the velocity 
field, we set V ang = and 

v rn (u = 0, r, x A ) = -(VTTE +^E + 2/r)(l - 2/r) (39) 

where v rn is V\ renormalized to be well behaved on the event horizon r = 2 (explicitly 
v rn = v i(l — 2/r) 2 ). E represents the energy at infinity of a unit mass particle, in the sense 
that at infinity If 1 ) 2 = E\ note that E > —1. In the tests described below we take 
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• R a = 4 and Rb = 7 

• A varying between 10 -5 and 10~ 13 

• G(x A ) = 1 (spherical shell), or G(x A ) given by Eq. fl38|) with /i = 0.4 (blob of matter 
with center at one of the poles whose density goes to zero at about 9 = 44°) 

• E = 2.25 

Analytically, the matter never reaches r = 2, but in practice, due to numerical diffusion, 
boundary conditions are still needed there; we impose 

P = 0, V ang = 0, v rn = 0. (40) 



B. Equation of state 

For dust, the equation of state is p = 0. We also investigate how the code copes with 
non zero pressure, and set up a test case with the density defined as above (with A = 10 -9 ), 
and with the equation of state 

p = 10- u p 1A . (41) 

Further, in order to keep p p (Eq. (0)) well-behaved for small p, we set p = whenever 
P < Pmaxl0~ 5 , where 

Pmax is the maximum value of p at u — 0. Although the magnitude of 
this pressure is rather small, it is encouraging to see that the code can handle it even though 
our implementation of the fluid equations is rather simple. More realistic equations of state 
would induce shocks that would recquire a sophisticated treatment of the matter equations. 
We defer this to future work. 

C. Physical interpretation of the initial data 

First, it is necessary to clarify the issue of units. We are using geometric units in which 
G = c = 1, and everything is given in terms of a unit of length. However, this unit is not 
metres, but rather it is the distance corresponding to a unit change in the radial coordinate 
r. We will call the geometric unit of length 1L, and in order to convert to S.I. units we need 
to know the value of 1L in metres. We write 

1L = ro [metres]. (42) 

For example, if we have a 10M Q black hole with event horizon at r = 2, then tq would 
be 14766 [metres]. We use the notation that quantities in geometric units will be denoted 
without suffix, and those in S.I. units will have the suffix 5. Then the conversions for speed 
v, mass M, density p and pressure p are: 
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v x 2.998 x 10 8 = v s [m/s] (43) 
M x 1.347r x 10 27 = M s [kg] (44) 

1.347 x 10 27 ri , „, 

P x ~2 = Ps [ k g/ m ] (45) 

'o 

1.210 x 10 44 

P x - 2 = Ps [Pa] (46) 

'o 

Note that 1.347 x 10 27 is c 2 /G, and 1.210 x 10 44 is c 4 /G in S.I. units. 

We now use these conversions to determine, in S.I. units, the various parameters describ- 
ing the initial data specified above, when the mass of the black hole is 1OM . For the case 
A = 10 -9 we find 

• Spherical shell: p s = 8.34 x 10 8 kg/m 3 , v = 0.90c, p s = or p s = 8.47 x 10 10 Pa, 
M = 3.37 x 1O" 6 M ; 

• Localized pulse: ps = 2.14 x 10 7 kg/m 3 , v = 0.90c, ps = or p s = 5.01 x 10 8 Pa, 
M = 6.11 x 1O" 9 M ; 

where v is the proper inward radial velocity, and the values of p$, ps and v are given at the 
point of maximum density at u = 0. For A ^ 10 -9 , ps and M scale linearly in A; ps scales 
as A 1 ' 4 ; and v is unchanged. 



D. Spherical collapse 

The first test consisted in the collapse of a spherical shell onto a black hole We set the 
black hole mass = 1, the inner shell R a = 4 and the outer shell Rb = 7. Since the resulting 
spacetime is spherically symmetric there should be no gravitational waves emitted by the 
system. We confirmed this behavior by monitoring the News function for different initial 
amplitudes A = 10~ (2n+5) (with n = 0..4) over time. The value of each polarization mode 
converges to zero as the discetrization gets refined. The matter collapses onto the black hole 
and the run proceeds smoothly. Figure [I] shows the evolution of p for the case with pressure. 

In this case it is fairly straightforward to calculate analytically the Bondi mass of the 
system, and using Maple we found that (for p = 0, A = 10~ 6 ) M = 1.0007395. By increasing 
the resolution, we checked that this value is approached, with second order accuracy, by the 
value computed with the code. For instance, with a grid of N x = 165 and = 65, the 
obtained value for the mass is M = 1.0007408; which agrees quite well with the expected 
value. 



E. Black hole-matter ball collision 

To study the collision of a dust ball-black hole system we again set the black hole mass 
M = 1 and the inner and outer radius of the blob with R a = 4 and R^ = 7; finally, the 
localization on the sphere is set by choosing /x = 0.4. This configuration describes a ball 
of dust with center at one of the poles that goes to zero at about 9 = 44° with compact 
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support in [4, 7] in the radial direction. Note that in this case, although we initially set 
Vang = 0, the self gravitational field of the dust ball induces a non zero V ang towards the 
pole. Figure displays the evolution of p for the case with pressure (A = 1CT 9 ). Again, the 
evolution proceeds smoothly as the "blob" of matter collapses on to the black hole and the 
plus polarization mode of the news is obtained at null infinity (see figure f|). 

Although an analytic solution to this problem is unknown, one can check consistency of 
the obtained simulation by observing that the cross polarization mode (A/" x ) must vanish 
(since the problem is still axysimmetric) . We confirmed this behavior, by increasing the 
number of grid points and plotting the logarithm of the norm of A/" x (at u = 0.15) vs. 
discretization size and observing its convergence to zero. As shown in figure the slope of 
1.9 is in agreement with second order convergence. The convergence test was performed with 
A = 1CT 9 and with non- zero pressure. We should add, however, that at later times (u > 10), 
the convergence rate is reduce to about 1.5. This is expected because of the non-centered 
scheme used at the black hole boundary (see section [TV]) . 

Another check made was that the path of the peak density should be a geodesic of the 
background spacetime. This was indeed found to be the case. For example, using Maple we 
found that at u = 2 r should be 4.9169; and numerically we found, for the case A = 10~ 7 
and non- zero pressure, that r = 4.9180 (for a grid with N x = 85 and iVg = 65). 

VI. CONCLUSION 

In this paper, we have incorporated matter, in the form of a perfect fluid, into the 
characteristic code for the Einstein field equations: we have found explicit forms for the 
various evolution and field equations, we have shown how these equations are discretized, 
and we have carried out a number of tests on the resulting code. The code is stable and 
convergent, and its validation has included the following tests 

• The peak of the matter "blob" follows a geodesic of the background spacetime. 

• The code is not written with any symmetries, yet symmetric initial data leads to the 
appropriate part of the gravitational radiation vanishing. 

• The initial mass of the system, as calculated by the code, agrees with the initial mass 
calculated analytically. 

This paper has not computed a solution to any real problem in astrophysics. Never- 
theless, we have shown that the present code is capable of modeling a variety of situations 
where matter is captured by a black hole; although our treatment of the hydrodynamics 
would need to be amended in order to be able to handle shock waves. Even so, the present 
code should be able to contribute towards an understanding of a number of problems in 
black hole astrophysics. 
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APPENDIX 



The evolution equations for the matter variables are given below. Note that V w is 1+W/r. 
V ang , u = ^({2e 2/3 (dp)V 2 ng + 2e 2 ^V ang (dp)V ang )K - AV ang p, r v r 2 - 2V ang p, r V ang Ur 2 

- AV ang p, uVl r 2 + 4e 2/3 (g P )r 2 - 2e 2l3 (Qp)V 2 ng J - 2V angVl (Qp)Ur 2 - 2V angVl (dp)Ur 2 
+ 4V ang p, rVl V w r 2 - 2V 2 ngP , r Ur 2 - 2e 2 f 5 V ang ($p)V ang j) 

+ ^r 2 ^v ir 2 UmV ang - e 2 ^dJ)V 2 ng - e 2 ^dJ)V 2 ng + 4^^(3/3)^ 

- AV ang>r v r 2 + W w V ang>rVl r 2 - 2e 2l) (W ang )V ang J - 2v ir 2 {W)V ang 

+ (2e 2 ?(W ang )V ang + 2e 2f3 (W ang )V ang )K - 4V w v 2 r 2 ((5[3) - 2e w (W ang )V ang J 
+ 2e 2 PjK(Q.J)V ang V ang + 2e 2fi JK(<5J)V ang V ang - 2 Vl U(W ang )r 2 

- 2e 2 \m)V ang V ang + 2vlr 2 (W w ) - Ae 2(S JJ(QK)V ang V ang - 2V ang , r UV ang r 2 

+ 8v v ir 2 (d(3) - 2 Vl U(dV ang )r 2 - 2r 2 Vl {W)V ang - 2V ang , r UV ang r 2 ) (47) 

eWjFM-Fgjp + p)) e 2 ^-F a + F oVlPp ) 

P ' U " v 2 M(P + P) ~ 1) ' ^ v 1 (p p (p + p) - 1) 1 " J 

with F and F a given by 
P + P 



Ae 2/3 r 



2 {[-2{m)V ang e 2 ^ - 2(QK)V ang e 2 ?) J J - 2(dJ)V ang e 2 ^ 



+ ((3J)W * + (<5J)V ang e 2p )KJ+ {$J)V ang e 2p + (<5J)V ang e 2p )KJ 
+ 8 Vl rV w - Av^.r 2 - 2 Vl r 2 {W) + (-2(W ang )e 2 ? - A(Q(3)V ang e 2l3 )J + AV^r 2 

- 2(dJ)V ang e 2(3 + (-A(B(3)V ang e^ - 2(W ang )e w )J 

+ {2{W ang )e 2 ? + 4(8/3)K„ 9 e 2 ^ + A(df3)V ang e 2 ^ + 2{W ang )e 2 ?)K 

- 2UV ang , r r 2 - 2UV ang , r r 2 + A Vl r 2 V w , r - 8v r - 2r 2 U, r V ang - 2U(Bv 1 )r 2 

- ArUV ang - ArUV ang - 2r 2 U r V ang - 2 Vl r 2 {W) - 2f/(5t; 1 )r 2 ) 

+ 4^2 ( - 2 P,rVan 9 Ur 2 - 2(Bp)V ang e 2 ^ J - Ap , r v r 2 - 2(dp)V ang e 2l3 J - 2p, r V ang Ur 2 

- 2 Vl (Qp)Ur 2 + (2(Qp)V ang e 2/3 + 2((5p)V ang e 2/3 )K + Ap^V w r 2 - 2v l (c5p)Ur 2 ) (49) 
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F - = iJ^s ( - ^ lP<r V ang Ur 3 - 2v 2 1 (dp)Ur 3 - 2v 2 1 (Bp)Ur 3 - 4v lP , r v r 3 + Av 2 lP>r V w r 3 

- 2rv l (3p)V ang e 2/3 J - 2rv 1 (dp)V ang e 2fS J + (2rv 1 (dp)V ang e 2fS + 2rv l (dp)V ang e 2 ^)K 

- 2v lPir V ang Ur 3 + 4p, r r 3 e 2/3 ) 

+ ^3 ((-2r(5^i)Kn 9 e 2/3 + 2V 2 ng e 2 ?) J + 8/3,^^ r 3 - 2^(6^ + 2v\r 3 V w , r 

- 2v 1)T UV ang r 3 + W an9 - 2 Vl r 3 U >r V ang + W w v^ r v x r 3 - 2v l)T UV ang r 3 
+ KJr J r Kn 9 Kn 9 e 2/3 - 2v 1 U(dv 1 )r 3 + Af3, rVl r 3 UV ang - Av hr v r 3 - rJ r V 2 ng e 2 ^ 
+ (-4V ang V ang e 2 P + 2r(g Vl )K n9 e 2/3 + 2r(B Vl )V ang e 2f3 )K - W w (3, r vjr 3 

- 2 Vl r 3 U r V ang + (-2r(d Vl )V ang e^ + 2V 2 ng e 2l3 )J - 2JJrK,V ang V ang e 2 ^ 

+ KJrJ r V ang V ang e 2 P - rJ r V 2 ng e 2 ^ . (50) 
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FIGURES 



u=0.0 u=6.3 




FIG. 1. Density profiles vs. time. The figure shows 4 different snapshots of p at the equator. 
The inner part "hole" in the pictures corresponds to the black hole r = 2m radius while the outer 
part corresponds to r = oo. The collapse of matter onto the black hole can be clearly seen. The 
plots are for the case A = 10~ 9 and p ^ 0. 
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A 



B 




FIG. 2. Density profiles vs. time for the localized "blob" of matter at times u = (A), u = 7.2 
(B), u = 14.3 (C) and u = 21.5 (D). The value of initial amplitude parameter for this case is 
A = 10~ 9 , and p ^ 0. The figure shows 4 different snapshots of p at y = (on the northern 
hemisphere) as a function of the compactified radial coordinate (x = 2/3 being the black hole 
radius and x = 1 corresponding to null infinity). The pulse collapses onto the black hole remaining 
localized. 
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-7.9 




-9.4 1 — -J 1 1 — — 1 1 1 

-2.2 -2.0 -1.8 -1.6 -1.4 

FIG. 3. Convergence of the cross polarization mode in an axysimmetric spacetime. The figure 
shows the logarithm of N x (u = 0.15) vs. the logarithm of the discretization size. The slope of 1.9 
is in good agreement with second order convergence. The plot is for the case A = 10~ 9 and p ^ 0. 
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FIG. 4. The "plus" component of the news on the northern hemisphere for the case A = 10 
and p / O.at times u = (A), u = 7.2 (B), u = 14.3 (C) and u = 21.5 (D). 
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